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Abstract 

We extend the technique of using the Trapezoidal Rule for efficient 
evaluation of the Special Functions of Mathematical Physics given by 
integral representations. This technique was recently used for Bessel 
functions, and here we treat Incomplete Gamma functions and the 
general Confluent Hypergeometric Function. 
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1 Introduction 



In a recent work I was able to demonstrate the power of the Trapezoidal Rule 
for numerical integration [1], 

/oo 
dxf(x)^Yl hf(nh), (1.1) 
-oo n 

applied to integral representations of the various Bessel functions. [2] Here, h 
is the mesh size, which will be systematically reduced to show the convergence 
of the method; and the sum over mesh points n, nominally from — oo to +oo, 
will be truncated when the contributions to the sum are below the desired 
accuracy. Now that general technique is applied to other members of the 
venerable family called the Special Functions of Mathematical Physics. 

First, we look at the Incomplete Gamma function (the earlier paper did 
deal briefly with the complete Gamma function); and we find that a partic- 
ular approach initiated some years ago by Talbot [3], employing the inverse 
Laplace transform, is very useful. 

Then, following more traditional lines, we start with the textbook integral 
representation for the general Confluent Hypergeometric function and change 
variables so that it behaves remarkably well under the Trapezoidal Rule. 

Lots of numerical results are presented, showing the rapid convergence 
and versatility of the calculational method. 



2 Incomplete Gamma Function 

The definiton is 

7(s,x) = / dtt 3 ^ 1 e _t ; (2.1) 
Jo 

and I want to follow Talbot's[3] idea, further developed by Weideman and 
Trefethen [4] [5], about looking at inverse Laplace Transforms. 
Taking the Laplace transform of (2.1), we find 

/ dxe-y 7 (s, x) = T(s) y- 1 (1 + y y s . (2.2) 
Jo 

The inverse transform gives us back the original function: 
l(s,x) = I ^Jdye x yy' 1 (l + y )- = £^ J ' dy y' 1 (1 + y/x)-, (2.3) 
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where we have assumed, for now, that X IS cl positive real number, so that 
the second equality above is just the result of scaling the integration variable 

y- 

The contour of this y-integral starts way out in the third quadrant and 
ends far out in the second quadrant, crossing the real axis at a positive value 
of y. This is the same contour we earlier used [2] for calculating the inverse 
of the complete Gamma function. We find that for x — > this gives us just 
that earlier formula; although the original definition, Eq. (2.1), gives the 
complete Gamma function when x — >■ oo. 

I have not seen this integral representation (2.3) for j(s,x) before. It 
should be great for using the Trapezoidal Rule. We start by choosing a 
simple contour: 

y = c + 1 — cosh(u) + i sinh(u), — oo < u < oo; (2.4) 

and c is where the contour crosses the real axis. We see that the integral 
(2.3) has a simple pole at y = and a branch point at y — —x. For c > we 
have the function 7(5, x); and if we take —x < c < 0, then we get the answer 
7(a,x)-r(a) = -XT dtr- 1 e'K 

Looking for the point of constant phase to be at this crossing, we find the 
formula 

c = [(s + 1 - x) + y/(s + 1 - x) 2 + Ax ]/2 (2.5) 

The data in Tables 1 and 2 below show some examples of the ratio P(s, x) = 
j(s,x)/r(s) computed by this method. These calculations, as in Ref[2], are 
done with standard double precision accuracy (10~ 16 ) and I truncate the sum 
over mesh points when the fractional contributions are less than this amount. 



Table 1. Computations of P(s,x) using the Trapezoidal Rule 



1/h 


P(0.1,l) x 10 


P(l,0.1) x 100 


P(0.1,0.1) x 10 


1 


9.85296 26362 19827 


9.58023 11961 90712 


8.37021 11048 74736 


2 


9.75973 88897 94659 


9.51192 39220 99238 


8.27666 75413 58269 


4 


9.75872 65150 00294 


9.51625 80387 73782 


8.27551 75852 50319 


8 


9.75872 65627 36719 


9.51625 81964 04048 


8.27551 75958 58505 


16 


9.75872 65627 36723 


9.51625 81964 04037 


8.27551 75958 58505 
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Table 2. Computations of P(s,x) using the Trapezoidal Rule 



1/h 


P(l,l) x 10 


P(10,10)xl0 


P(1000, 1000) x 10 


1 


6.35418 14003 89701 


5.69571 39376 33220 


5.45885 48876 16142 


2 


6.32017 05027 55139 


5.41549 12109 30272 


5.11064 21436 17747 


4 


6.32120 56442 73888 


5.42070 15780 49574 


5.04658 65045 16275 


8 


6.32120 55882 85574 


5.42070 28552 84053 


5.04204 40307 32861 


16 


6.32120 55882 85577 


5.42070 28552 81477 


5.04205 21812 85238 


32 




5.42070 28552 81479 


5.04205 24418 02230 


64 






5.04205 24418 02222 



The number of mesh points used for the last line of data in Table 1 was 
153,151,153 respectively; and for Table 2 the numbers were 151, 285, and 
841. 

Special cases of the Incomplete Gamma Function include the Error Func- 
tion and various Exponential Integrals. 



3 Confluent Hypergeometric Functions 

The focus here is on the integral representation, 

C(a,b;x) = Cdtf 1 - 1 (1-tf- 1 e xt , a,b>0. (3.1) 
J o 

This is connected to the traditional definition of the Confluent Hypergeo- 
metric Function, as follows. 

iFi(a, c = b + a;x) = M(a, c = b + a;x) = B(a, by 1 C(a, b; x). (3.2) 

And this introduces the classical Beta-function, 

B(a, b) = T(a) r(6)/r(a + b) = C(a, b;x = 0). (3.3) 

where that last identity comes from the fact that the Hypergeometric func- 
tions are defined to have the value 1 at x — 0. 

We propose to evaluate these functions by the Trapezoidal Rule, after 
making some convenient (real) coordinate transformations under the integral: 

t = hl + tanh(u)) = e u /( y e u + e- u ); u = sinh(v); (3.4) 
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and then we apply the Trapezoidal rule to the infinite integral over v. 
Some numerical results are shown in Tables 3-7 below. 



Table 3. Computations of C(a — 1., b — 1.; x) 



1/h 


(x = 0.) 


(x = 1.) 


(x = 100.) x 10" 41 


1 


1.00107 86347 90329 


1.72266 65011 24113 


1.34247 60795 12472 


2 


1.00000 01397 11616 


1.71828 29887 33384 


2.83640 72153 82483 


4 


1.00000 00000 00001 


1.71828 18284 59079 


2.68831 33551 31161 


8 




1.71828 18284 59044 


2.68811 71418 07843 


16 






2.68811 71418 16129 


32 






2.68811 71418 16129 



Table 4. Computations of C(a = 0.1, b — l.;x) 



1/h 


(x = 0.) 


(x=l.) 


(x = 100.) x 10" 41 


1 


9.99991 08876 77476 


11.21464 29773 5867 


1.34425 85914 11300 


2 


10.00000 00993 10149 


11.21300 57977 0100 


2.86472 18192 35488 


4 


9.99999 99999 99998 


11.21300 52032 3319 


2.71295 94923 56232 


8 


9.99999 99999 99998 


11.21300 52032 3318 


2.71278 37414 62587 


16 






2.71278 37414 71210 


32 






2.71278 37414 71210 



Table 5. Computations of C(a — 0.1, b — 10.; x) 



1/h 


(a: = 0.) 


(x = 1.) 


(a; = 100.) x 10~ 29 


1 


7.53951 55733 97154 


7.63169 79919 69269 


3.40378 23912 96282 


2 


7.59131 58970 84419 


7.67046 16609 55571 


1.70396 29450 35887 


4 


7.59138 00009 01282 


7.67049 54154 30269 


1.09112 36709 31219 


8 


7.59138 00009 11013 


7.67049 54154 32864 


1.07365 08998 41708 


16 


7.59138 00009 11017 


7.67049 54154 32878 


1.07365 07978 79342 


32 






1.07365 07978 79343 
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Table 6. Computations of C(a = 10., b = 0.1; x ) 



1/h 


(x = 0.) 


(x = 1.) 


(x = 100.) x 10~ 44 


1 


7.53951 55733 97154 


20.26525 88483 9334 


1.69844 72080 59734 


2 


7.59131 58970 84419 


20.44048 93266 6513 


1.59374 46727 15385 


4 


7.59138 00009 01285 


20.44076 89724 1024 


1.59966 19996 26905 


8 


7.59138 00009 11014 


20.44076 89724 7923 


1.59966 08127 76379 


16 


7.59138 00009 11021 


20.44076 89724 7924 


1.59966 08127 76238 


32 






1.59966 08127 76246 



Table 7. Computations of C(a = 0.1, b = 0.1; x) 



1/h 


(x = 0.) 


(x = 1.) 


(x = 100.) x 10~ 44 


1 


19.71352 25754 7283 


35.95446 58322 5812 


1.70487 63130 76804 


2 


19.71463 97119 7631 


35.95643 50355 3144 


1.61024 26167 51076 


4 


19.71463 94890 5016 


35.95643 47587 2009 


1.61504 43786 53350 


8 


19.71463 94890 5015 


35.95643 47587 2014 


1.61504 16242 89936 


16 




35.95643 47587 2013 


1.61504 16242 89858 


32 






1.61504 16242 89859 



The number of mesh points needed to obtain the best results shown in 
the Tables above was always under 200. 

Looking at this data, one is reminded of a particular virtue of the present 
method. [1] The integral in Eq. (3.1) has, for non-integral values of a and b, 
singularities at the end points t — 0,1. This means that traditional methods 
of numerical integration (Simpson's rule, Richardson extrapolation, Gauss 
quadrature, etc.) would not work at all well. But the transformation Eq. 
(3.4) used for our approach with the Trapezoidal Rule, handles those sin- 
gularities very nicely: they are moved to infinity and smothered beneath 
exponential decays. 

The same method used above should be applicable to the general Hyper- 
geometric function, 

2 F 1 (a, b; c; z) = Bib, c - by 1 C dt t^ 1 (1 - t) c ~ b - 1 (1 - zt)~ a , (3.5) 

Jo 

although I have not done those calculations. There are various transforma- 
tion formulas for 2 F\ that let one deal with the singularities when z ap- 
proaches 1 or infinity. 
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4 Discussion 



The numerical results displayed above show very nice rates of convergence 
using the Trapezoidal Rule. The programming is straightforward; and there 
is freedom for the user to explore various avenues. 

While some other authors, e.g. [4], have focused on finding optimal con- 
tours, my own opinion is that this general method is so powerful and robust 
that such refinements may be more a distraction than a benefit. 

Since I have restricted myself here to real values of the coordinates (x) 
and the parameters (s,a,b), it is left for others to explore the extension of 
those variables into the complex planes. 

I am grateful to J.A.C.Weideman for some stimulating conversations. 
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